#include "structMeshFunctions.hpp"
#include "rectMeshFunctions.hpp"
#include "boundaryPatch.hpp"
#include "uniMesh.hpp"
#include "rectMesh.hpp"
#include "structMesh.hpp"

#include <stdexcept>


namespace hms {


/**
 * @brief Returns the boundary location of an edge according to its index and
 * the provided mesh size.
 * 
 * @param edge The edge index.
 * @param meshSize The number of cells in topological x- and y-direction.
 * @return BoundaryLocation Can be bottom, top, left or right.
 */
BoundaryLocation whereIs( Index edge, const Ref<const Array2i>& meshSize ){
	if ( edge < meshSize.x() )
		return BoundaryLocation::bottom;
	else if (
		edge >= meshSize.prod() &&
		edge <  meshSize.x() * (meshSize.y()+1)
	)
		return BoundaryLocation::top;
	else {
		edge -= meshSize.x() * (meshSize.y()+1);

		if ( edge % (meshSize.x()+1) == 0 )
			return BoundaryLocation::left;
		else {
			assert( (edge+1) % (meshSize.x()+1) == 0 &&
			"Edge is not at boundary!" );
			return BoundaryLocation::right;
		}
	}
}

/**
 * @overload
 * 
 * Uses the boundary patch's first edge index to determine the location.
 * Boundary patch must not be empty.
 */
BoundaryLocation whereIs(
	const BoundaryPatch& patch,
	const Ref<const Array2i>& meshSize
){
	assert( !patch.empty() && "Boundary patch is empty!" );
	return whereIs( patch.edges()(0), meshSize );
}


Array<BoundaryLocation, 1, Dynamic> whereAre( 
	const Ref<const ArrayXi>& edges,
	const Ref<const Array2i>& meshSize
){
	Array<BoundaryLocation, 1, Dynamic> locations (1, edges.size());

	for (int i = 0; i< edges.size(); i++) {
		locations(i) = whereIs(edges(i), meshSize);
	}

	return locations;
}


BoundaryLocation whereIsVertex( Index vert, const Array2i& meshSize){
	using BL = BoundaryLocation;

	const Index
		& nx { meshSize.x() },
		& ny { meshSize.y() };

	if ( vert < nx + 1 )
		return BL::bottom;
	else if (
		vert >= (nx+1) *  ny    &&
		vert <  (nx+1) * (ny+1)
	)
		return BL::top;
	else if (  vert    % (nx+1) == 0 )
		return BL::left;
	else if ( (vert+1) % (nx+1) == 0 )
		return BL::right;
	else
		throw std::runtime_error( "Starting vertex not on boundary!" );
}


Index nHorizontalEdges( Index nx, Index ny ){
	return nx * (ny + 1);
}
Index nHorizontalEdges( const Ref<const Array2i>& meshSize ){
	return nHorizontalEdges( meshSize.x(), meshSize.y() );
}


/**
 * @brief Returns true if the edge index resolves to a horizontal edge according
 * to the provided mesh size.
 * 
 * @param edgeIdx Edge index within the provided mesh size.
 * @param meshSize Mesh size according to which it is determined, whether or
 * not the edge index refers to a horizontal edge or not.
 */
bool isHorizontal( Index idx, const Ref<const Array2i>& meshSize ){
	return idx < nHorizontalEdges(meshSize);
}


template<typename MeshType>
auto whichRange( const StructMeshBase<MeshType>& mesh, Index edge ) ->
	std::pair<EdgeRange, Index>
{
	assert( 0 <= edge && edge < mesh.nEdges() );
	if ( edge < mesh.nCellsX() * (mesh.nCellsY()+1) )
		return {EdgeRange::horizontal, edge};
	else
		return {EdgeRange::vertical, edge - mesh.nCellsX()*(mesh.nCellsY()+1)};
}


template<typename MeshType>
Array2i parentCells(
	const StructMeshBase<MeshType>& mesh,
	Index edge,
	Index order
){
	assert( edge >= 0 && edge < mesh.nEdges() && "Edge index out of bounds!" );
	const Index
		nx { mesh.nCellsX() },
		ny { mesh.nCellsY() },
		nh { nHorizontalEdges(nx, ny) };

	static constexpr Index ng {1};
	
	Array2i sides { -order, order-1};

	if (edge < nh)
		return (nx + 2*ng) * (sides + ng + edge / nx) + ng + edge % nx;

	edge -= nh; // now vertical edges start at zero again
	return sides + (nx + 2*ng) * ( ng + edge / (nx+1) ) + ng + edge % (nx+1);
}

template<typename MeshType>
void addStandardBoundaryPatches( StructMeshBase<MeshType>& mesh ){

	mesh.addBoundaryPatch( "bottom",
		ArrayXi::LinSpaced( mesh.nCellsX(), 0, mesh.nCellsX()-1 ) );
		
	mesh.addBoundaryPatch( "top",
		ArrayXi::LinSpaced( mesh.nCellsX(),
			mesh.nCellsX()*mesh.nCellsY(),
			mesh.nCellsX()*(mesh.nCellsY()+1) - 1) );

	mesh.addBoundaryPatch( "left",
		ArrayXi::LinSpaced( mesh.nCellsY(),
			mesh.nCellsX()*(mesh.nCellsY()+1),
			2*mesh.nCellsX()*mesh.nCellsY() + mesh.nCellsY() - 1 ) );
			
	mesh.addBoundaryPatch( "right",
		ArrayXi::LinSpaced( mesh.nCellsY(),
			mesh.nCellsX()*(mesh.nCellsY()+2),
			mesh.nEdges()-1 ) );
}



Array4i structMeshCorners( const Ref<const Array2i>& meshSize ){
	Index
		nx { meshSize.x() },
		ny { meshSize.y() };

	return {
		0,              // bottom left
		nx,             // bottom right
		(nx+1)*ny,      // top left
		(nx+1)*(ny+1)-1 // top right
	};
}


Index nextBoundaryVertex(
	Index beg,
	bool clockwise,
	const Array2i& meshSize
){
	using BL = BoundaryLocation;
	using C = Corner;

	Array4i corners { structMeshCorners(meshSize) };
	C corner {C::nCorners};
	for ( Index i{0}; i<4; ++i ){
		if ( beg == corners[i] ){
			corner = static_cast<C>(i);
			break;
		}
	}

	BL loc { whereIsVertex( beg, meshSize ) };

	assert( (
		   ( corner == C::bottomLeft  && (loc == BL::bottom || loc == BL::left ) )
		|| ( corner == C::bottomRight && (loc == BL::bottom || loc == BL::right) )
		|| ( corner == C::topLeft     && (loc == BL::top    || loc == BL::left ) )
		|| ( corner == C::topRight    && (loc == BL::top    || loc == BL::right) )
		|| corner == C::nCorners
	) );

	Index next {signalCode::uninitialised};

	if (
		(!clockwise && loc == BL::bottom && corner != C::bottomRight) ||
		( clockwise && loc == BL::top    && corner != C::topRight   )
	)
		next = beg + 1;
	else if (
		(!clockwise && loc == BL::top    && corner != C::topLeft    ) ||
		( clockwise && loc == BL::bottom && corner != C::bottomLeft )
	)
		next = beg - 1;
	else if (
		(!clockwise && (loc == BL::right || corner == C::bottomRight)) ||
		( clockwise && (loc == BL::left  || corner == C::bottomLeft ))
	)
		next = beg + (meshSize.x() + 1);
	else if (
		(!clockwise && (loc == BL::left  || corner == C::topLeft )) ||
		( clockwise && (loc == BL::right || corner == C::topRight))
	)
		next = beg - (meshSize.x() + 1);
	else
		throw std::logic_error(
			"Implementation issue: branch should be unreachable."
		);

	return next;
}



template<typename MeshType>
bool validCellIndex( const StructMeshBase<MeshType>& mesh, Index i, Index j ){
	/* i-indices can go from -nGhostLayers to nCellsX-1+nGhostLayers, and
	 * analogous for j, however, the corner ghost cells (e.g. i=-1 and j=-1)
	 * don't exist. */
	Index nx { mesh.nCellsX() };
	Index ny { mesh.nCellsY() };
	static constexpr Index ng {1};
	return (
		(
			i >= -ng && i < nx+ng  && /* full range for i, but j constrained */
			j >=  0  && j < ny
		) || (
			i >=  0  && i < nx     && /* or the other way around */
			j >= -ng && j < ny+ng
		)
	);
}


template<typename MeshType>
bool validCellIndex( const StructMeshBase<MeshType>& mesh, Index cell ){
	/* mustn't be negative and mustn't exceed number of cells.
	 * Furthermore, must not be one of the corners, which only border ghost
	 * cells. Non-negative and non-zero (corner) means greater than zero, 
	 * smaller than number of cells and not a corner means smaller than
	 * (m+2)*(n+2)-1. */
	Index nx  { mesh.nCellsX() };
	Index ny  { mesh.nCellsY() };
	Index ng { mesh.nGhostLayers() };
	return
		/* within bounding box */
		   cell > 0
		&& cell < (nx+2*ng)*(ny+2*ng)
		/* when in left or right ghost region, then not in top or bottom one and
		 * when in top or bottom ghost region, then not in left or right one */
		&& (
			( cell % (nx+2*ng) >= ng && cell % (nx+2*ng) < nx + ng ) ||
			( cell / (nx+2*ng) >= ng && cell / (nx+2*ng) < ny + ng )
		);
}


/* functions that need to be redefined for release mode */
#ifndef NDEBUG

template<typename MeshType>
void assertValidCellIndex( const StructMeshBase<MeshType>& mesh, Index cell ){
	assert( validCellIndex(mesh, cell) && "Cell index invalid!" );
}

template<typename MeshType>
void assertValidCellIndex(
	const StructMeshBase<MeshType>& mesh, Index i, Index j
){
	assert( validCellIndex(mesh, i, j) && "Cell index invalid!" );
}

#else // NDEBUG

/* these functions shall be optimised away in release mode, otherwise
 * the compiler sees unused variables. */
template<typename MeshType>
void assertValidCellIndex( const StructMeshBase<MeshType>&, Index ){}

template<typename MeshType>
void assertValidCellIndex( const StructMeshBase<MeshType>&, Index, Index ){}

#endif


#define INSTANTIATE_FOR_STRUCTURED_MESHTYPES(MeshType) \
template Array2i parentCells( const StructMeshBase<MeshType>&, Index, Index ); \
template std::pair<EdgeRange, Index> whichRange( const StructMeshBase<MeshType>&, Index ); \
template void addStandardBoundaryPatches( StructMeshBase<MeshType>& ); \
template bool validCellIndex( const StructMeshBase<MeshType>&, Index ); \
template bool validCellIndex( const StructMeshBase<MeshType>&, Index, Index ); \
template void assertValidCellIndex( const StructMeshBase<MeshType>&, Index ); \
template void assertValidCellIndex( const StructMeshBase<MeshType>&, Index, Index );

#include "structMeshInstantiator.hpp"

} // namespace hms